Biological materials.

Adult Lamellibrachia luymesi specimens were collected from seep localities in the Mississippi Canyon at 754 m depth in Gulf of Mexico (N 28°11.58’, W 89°47.94’), using the R/V Seward Johnson and Johnson Sea Link in October 2009. All samples were frozen at 80˚C following collection.

Genome sequencing and assembly.

Vestimentum tissue was dissected from one individual of worm, and high molecular weight genomic DNA was extracted using the DNeasy Blood & Tissue Kit (Qiagen) according to the manufacturer’s protocols. Sequencing a total of six paired-end or mate-pair genomic DNA libraries with insert sizes ranging from 180 bp to 7 kb were performed by The Genomic Services Lab at the Hudson Alpha Institute in Huntsville, Alabama on an Illumina HiSeq 2000 platform (see details in Table S1). Paired-end libraries (180 bp, 400 bp, 750 bp) were prepared using the 125 bp TrueSeq protocols, and mate-pair libraries (3-5 kbp, 5-7 kbp) were generated using the Illumina Nextera Mate Pair Library Kit followed by size selection. Moreover, a 10X sequencing library was constructed using the 10X Chromium protocol (10X genomics) at the Hudson Alpha Institute. The finished library was sequenced on an Illumina HiSeqX platform, using paired 151 bp reads with a single 8 bp index read.

Our workflow of genome assembly was shown in Fig. S1. Two The paired-end and 10X raw reads were checked using FastQC v0.11.5 (Andrews and others 2010) and quality filtered (Q score >30) using Trimmomatic v0.36 (Bolger, Lohse, and Usadel 2014). The estimation of genome size, level of heterozygosity and repeat contes of the Lamellibrachia genome was determined by analaysing the kmer histograms generated from the paired-end libraries using Jellyfish v2.2.3 (Marçais and Kingsford 2011) and GenomeScope (Vurture et al. 2017) (Fig. S2). The Mate-pair reads were trimmed and sorted using NxTrim v0.3.1 (O’Connell et al. 2015) which can recgonize and trim the artificial Nextera mate-pair circulation adapters. Only reads from category “mp” (true mate-pair reads) and “unkonwn” (mostly large insert size reads) were used for downstream scaffolding anlaysis. Reads from “pe” (paired-end reads) and “se” (single ends) categories were discarded.

Given that high heterozygosity of Lamellibrachia genome, all reads were assembled using Platanus v1.2.4 (Kajitani et al. 2014) with a kmer size of 32. Scaffolding was conducted by mapping Illumina paired-end and mate-pair reads to contigs genrated by Platanus using SSPACE v3.0 (Boetzer and Pirovano 2014). Gaps in the scaffolds were then filled with GapCloser v1.12 (Luo et al. 2012). Redundant allele scaffods were further remvoed using Redundans v0.13c with default settings (Pryszcz and Gabaldón 2016). Genome assembly quality was assessed using QUAST v3.1 (Gurevich et al. 2013). Completeness of obtained genome was assessed using BUSCO v3(Waterhouse et al. 2017) with Metazoa_odb9 database (978 busco genes).

Transcripome assembly and analysis

Total RNA was extracted from the plume, vestimentum and trophosome tissue from the same indivdisual of Lamellibrachia spicemen using Trizol. RNA-seq of adult tissues from plume, vestimentum and trophosome was performed using Illumina HiSeq 2000 platform in Hudson Alpha. After quality checking and trimming of raw sequencing reads, transcripts were assembled de novo with Trinity. Transcript isoforms with high similarity (≥ 95%) were removed with CD-HIT-EST v4.7 (Li and Godzik 2006). Transcript abundance was estimated with Bowtie v2.2.9 (Langmead and Salzberg 2012) and RSEM v1.2.26 (Li and Dewey 2011) by mapping reads back to the transcriptomic assembly based on transcripts per million (TPM). A tissue specifically expressed gene was defined as a gene that had over 75% of the total transcripts in a particular tissue based on TPM (Albertin et al. 2015). GO enrichment analysis was performed based on the GO annotation using PANNZER2 (Törönen, Medlar, and Holm 2018). Statistically overrepresented GO terms of trophosome-specific genes were identifed using Lamellibrachia gene models as background with AgriGO (Du et al. 2010).

Genome annotation.

Our genome annotation workflow was shown in Fig. S3. Gene models of Lamellibrachia genome were constructed following the Funannotate pipeline 1.3.0 (https://github.com/nextgenusfs/funannotate). Briefly, repeptive regions in the Lamellibrachia genome were identified usning RepeatModeler v1.0.8 (Smit and Hubley 2008) and were subsequently soft-masked using RepeatMasker v4.0.6 (Chen 2004). RNA-Seq data from different tissue were leveraged to improve the accuracy of gene prediction. RNA-Seq data were assembled de novo into transcriptomes using Trinity v2.4.0 (Haas et al. 2013) and HISAT 2.1.0 (Kim, Langmead, and Salzberg 2015) was used to algin RNA-Seq reads to the Lamellibrachia assembly. Transcrptome assemblies were then passed to PASA pipeline v2.3.3 (Haas et al. 2003) to identify high quality gene models. The aligned RNA-Seq data wes used to train the ab initio gene predictions using AUGUSTUS v3.3 (Stanke et al. 2006). Protein alignements from the SwissProt database to “Lamellibrachia” assembly were generated using exonerate (Slater and Birney 2005) and Trinity/PASA transcripts were aligned to the genome using Minimap2 v2.1 (Li 2018). The tRNA genes were identified using tRNAscan-SE v1.3.1 (Lowe and Eddy 1997). Finally, EvidenceModeler 1.1.0 (Haas et al. 2008) was used to combine all the evidences of gene prediction from protein alignemnts, transcritp alignments, and ab initio predictions to construct high quality gene models. Finally, functional annotations of predicted gene models were analyzed using several curated databases. KEGG orthology was assinged using the KEGG Automatic Annotation server. Gene models were further annotated with domain structure and protein identity by InterProScan (Zdobnov and Apweiler 2001) and SwissProt database, respectively. Secreted proteins were predicted using SignalP (Petersen et al. 2011) and Phobius (Käll, Krogh, and Sonnhammer 2007) using InterProScan.

Phylogenomics adn Molecular clock analysis

Analysis of siboglinid phylogeny was conducted from publically available 16 transcriptomic data in conjuction with Lamellibrachia proteome. Sequence assembly, annotation, homology evaluation, gene tree construction, parsing of genes trees to OGs, and supermatrix construction were conducted with Agalma (Dunn, Howison, and Zapata 2013). The final analyses presented here contained 13 siboglinids and four outgroups based on current understanding of annelid phylogeny. Maximum likelihood analyses were performed in IQTree under the best-fitting models for associated partition schemes determined by Modelfinder implemented in IQTree 1.6.3 (Nguyen et al. 2014) with ultrafast bootstrapping of 1000 replicates.

For molecular clock analysis, a relaxed molecular clock with a lognormal distribution and a Yule tree model was used in BEAST 2 v2.5.1 (Bouckaert et al. 2014). A calibration was placed on the node representing the most recent common ancestor (MRCA) of Osedax using a normal distribution with a mean of 100 MYA and a standard deviation of 10 following the findings of (Danise and Higgs 2015).Another calibrartion was placed on the node of MRCA of Serpulida and Sabellida using a normal distribution with a mean of 267 MYA (Sanfilippo et al. 2017). Molecular clock analyses with BEAST 2 consisted of two independent runs with 1 million MCMC generations sampled every 1000 generations. Convergence was checked and confirmed by comparing trace plots in Tracer making sure the effective sample size of each parameter was greater than 100 and that stationarity appeared to have been achieved. Log and tree fiels were combined using Logcombiner. A maximum clade credibility tree with mean heights was calculated using TreeAnnotater. The resulted time-calibrated tree was plotted using R package, phyloch, strap and OutbreakTools. Bayesian inference using a molecular clock resulted in identical branching patterns as analysis with IQTree.

Gene family analysis

After all-to-all Diamond v1.0 (Buchfink, Xie, and Huson 2014) BLASTP searches against 12 selecte lophotrochozoan proteomes, orthology groups (OGs) were identified using Orthofinder with a dfault inflation parameter (I=1.5). Gene ontology annotation was performed using PANTHER v13.1 (Mi et al. 2016) with the PANTHER HMM scoring tool (pantherScore2.pl). Gene family expansion and contraction was estimated using CAFÉ v2.1 (De Bie et al. 2006). For each gene family, CAFÉ generated a family-wide P value, with a significant P value indicating a possible gene-family expansion or contraction event. A branch-specific P value was also generated for each branch/node using the Viterbi method. A family-wide P value less than 0.01 and a branch/node Viterbi P value less than 0.001 was considered as a signature of gene family expansion/contraction for a specific gene family and specific species, respectively, as suggested in manual.

Hemoglobin and immunity-related genes.

Genome features

Results from high-throughput sequencing and genome assembly for Lamellibrachia are presented in Table # with at least 300 fold coverage using a combination of Illumina paired-end, mate-pair and 10X genomic sequencing with a variety of insert sizes. The haploid genome assembly size is ~687 Mb, similar to the result of genome size estimation using GenomeScope (~646 Mb). The N50 value of the assembled scaffolds and contigs is 373 Kb and 24 Kb, respectively. Although N50 lengths and assembly quality of Lamellibrachia are comparable to those of other annelids (e.g. Capitella teleta, Hebdella robusta), the overall genome completeness measured by BUSCO (~ 95%) is one of the highest among other lophotrochozoans (Table #), indicating the completeness of the genome assembly. With the support of RNA-seq data from different tissue, we estimated Lamellibrachia genome contains 38,998 gene models. The genome also exhibit similar level of heterozygosity (0.6%) and repeptitive seuqences (36.92%), compared to other lophotrochzoan genomes (Table #). Comparison among three annelid genome revealed

Hemoglobin evolution

Remarkably siboglinid Hbs are able to bind oxygen and sulfide simultaneously and reversibly at two different sites. The binding characteristics of these Hbs not only protect against the inhibition of aerobic respiration but also facilitate the delivery of sulfide and oxygen to the chemoautotrophic endosymbionts (Arp and Childress 1981). Siboglinids possesses a multihemoglobin system with three different extracellular hemoglobins (Hbs; V1, V2, and C1): two dissolved in the vascular blood, V1 and V2, and one in the coelomic fluid, C1 (Arp and Childress 1981, @zal1996multi). Siboglinid Hbs consist of four heme-containing chains (A1, A2, B1, B2) and four linker chains (L1-L4). Sulfur-binding capabilities are hypothesized to be due cysteine residues at key positions in Hbs, especially in the A2 and B2 chains. V1 Hb can form persulfide groups on its linker chains, a mechanism that can account for the higher sulfide-binding potential of this Hb (Zal et al. 1997). Moreover, part of the H2S binding affinity has also been hypothesized to be mediated by the zinc moieties bound to amino acid residues at the interface between pairs of A2 chains (Flores et al. 2005).

Phylogeny and molecular clock of tubeworm

Two additional Osedax taxa were added in the study compared to the previous siboglinid phylogenomic analyses (Li et al. 2017) (Table #). The final supermatrix dataset contains 191 genes single-copy orthologs. Bayesian inference with a relaxed molecular clock (Fig. 2) recovered the same topology as ML analysis in IQTree with strong nodal support (Fig. #). Both analyses strongly support Osedax is most closely related to the Vestimentifera + Sclerolinum clade and Frenulata is the early diverging group, as recently reported (Li et al. 2015). Within Vestimentifera, Lamellibrachia is sister to the remaining sampled vestimentiferans.

Molecular clock analyses based on phylogenomic dataset suggest modern siboglinid diversity originated in Mesozoic (223MYA ± 80 MY), conflicting with previous hypotheses indicate a Late Mesozoic or Cenozoic approximately 50-126 MYA. The previous analyses are solely based on estimation of nucleotide divergence of COI sequences on limited taxa sampling (mainly vestimentiferans). However, mitochondrial genes of vestiemntiferans may have experienced a “slow-down” in rate of nucleotide subsitution relative to other siboglinid lineages (???). Moreover, molecular clock analysis suggest a yonger split of Vestimentifernas during the Cenozoic (60 MYA). Recent analyses suggested that Jurassic tubes from Figueroa deposits is likely to have been made by vestimentiferans (Georgieva et al. 2017), but this hypothesisis not supported in our analysis. Thus, we provide the molecular clock of the siboglinid phylogenetic tree, placing a common siboglinid ancestor as far back as the Jurassic. Siboglinid vestimentifernas might exploit and adapt to vents and seeps during Cenozoic (<65 MYA). This adds to the growing evidence that the Cenozoic was a key period for the radiation of most dominant invertebrate taxa now occupying in deep-sea chomosynthetic communities (Vrijenhoek 2013).